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Abstract 

We examine the Melnikov criterion for transition to chaos in case of a single degree 
of freedom nonlinear oscillator with the Ueda well potential and an external periodic 
excitation term. Using effective Hamiltonian we have examined homoclinic orbits 
and cross-sections of stable and unstable manifolds which gave the condition of 
transition to chaos through a homoclinic bifurcation. 
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1 Introduction 



The Melnikov method has become a classical approach for predicting chaotic 
bechaviour in presence of saddle points on the basis of cross-sections of stable 
and unstable manifolds [1,2,3]. Usually, it is applied explicitly to systems which 
possess homoclinic orbits in multiple well potential like Duffing double well, or 
pendulum systems [2,3,4], or to the single well systems with a smooth poten- 
tial barrier against an unstable solution [5,6,7]. By signaling global homoclinic 
transition it builds a condition for creation of fractal boundary between at- 
traction basins and chaos appearance provided that a vibration amplitude is 
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large enough to reach this boundary. This is the reason why this condition can 
be easiest fulfilled in the region of nonlinear resonance. Comparing to other 
approaches [8,9,10], the above scenario is so clear and instructive that the 
main stream of research has been focused on smooth multi- well potentials 
where the basins of attractions belong to separate wells. 

In the present paper we will adopt the Melnikov method to an effective sys- 
tem described by double solutions. We will also investigate a possible fractal 
smearing of the basins of their attractions. Identifying a saddle point we will 
find homoclinic orbits there and finally define the corresponding Melnikov 
criterion. 
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Fig. 1. Different possible potentials of a single degree freedom system 
V(x) = V s i(x), V S 2(x), Vdi(x) and Vd2(x) as defined in Eq. 2 for 7 = 1.0 and 
<5 = 1.0. 

We start our analysis with one of the best known examples exhibiting chaotic 
solution, namely, the Ueda single well system 



where x is displacement ax is linear damping, \i sin Vtt is an external excitation 
while — 7X 3 is a cubic restore force (7 > 0). 

The above example is known from the pioneering Ueda work on chaotic sys- 
tems [11]. Note that there are substantial differences between a single Ueda 
well potential V^i without a linear term and a double well Duffing potential 
Vdi or an upside-down reflected double well potential Vd2 = —Vdi- On the 
other hand V s \ resembles Duffing potential with hard stiffness V S 2 
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Fig. 2. Phase portraits (Fig. 2a, c - with lines) and Poincare maps (Fig. 2a, c - with 
points) and corresponding time histories x{t) (Fig. 2b, d) for two sets of system (Eq. 

1) parameters: Fig. 2a,b with 7 = 1.0, /i = 7.5, a = 0.05, O = 1.0, Fig. 2c, d with 
7 = 1// = 0.1, a = 0.05, f2 = 0.21. Top Lyapunov exponents for cases Fig. 2a, b and 
Fig. 2c, d are Ai = 0.111 and Ai = 0.021, respectively. The initial conditions used in 
both cases x(t = 0) = 0.0, v{t = 0) = 2.5. 

where 7 and 5 are denned positive. In case of Vdiix) and Vd2(x) potentials have 
well denned inflexion points and extrema which correspond to unstable saddle 
fixed points (Fig. 1) while potentials V s i(x) and V s %{x) are fundamentally dif- 
ferent with a single minimum at x = and without any inflexion points. Note 
that previous applications of Melnikov theory based on existence of multiple 
extrema of type Vdi(x) or Vd2{x). 

On the other hand, Chakraborty, in his recent paper [12], noticed that one can 
expect to apply the Melnikov theory even for a single well potential V S 2 (Eq. 

2) after defining a new coordinate system. Motivated by this conjecture we 
are going to construct the Melnikov function and derive a necessary condition 
to system transition into a chaotic motion region. 

To explore this possibility further we have performed numerical simulations 
to find the regions of chaotic solutions (Eq. 1). In Fig. 2 we show the phase 
portraits and Poincare maps as well as time histories of chaotic solutions. Note 
Figs. 2a-b correspond to the original Ueda system [11], where 7 = 1.0, \i = 7.5, 
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Fig. 3. Top Lyapunov exponents versus f2 G [0.1 : 1.9] for two sets of system 
parameters; 7 = 1.0, fi = 7.5, a = 0.05 (Fig 3a) and Fig. with 7 = 1.0 /i = 0.1, 
a = 0.05 (Fig. 3b). For the smallest O the initial condition we assumed to be 
x(t = 0) = 0.0, v(t = 0) = 2.5 while for any next larger Q the final state values of 
x and v have been used as the initial conditions. 



a = 0.05, f2 = 1.0. Unfortunately, the large value fi makes any perturbation 
method non-relevant while Note Figs. 2c-d, show the chaotic solution for other 
choice of system parameters: 7 = 1.0, fi = 0.1, a = 0.05, Q = 0.21. These 
chaotic solutions are characterized not only by the fractal strange attractors 
and non-periodic time series (Figs. 2a-d) but also by positive top Lyapunov 
exponents Ai = 0.111 and Ai = 0.021 (for /i = 7.5 and fi = 0.1 respectively). 
In this second case (2c-d) a perturbation expansion in terms of fi could be 
performed. In spite of different system parameters the attractors (given by 
Poincare maps Figs. 2a and 2c) look similar. To show the regions of chaotic 
solutions we have plotted the top Lyapunov exponent as a function of Q in 
Fig. 3. For larger fi we observe the relatively wide region of chaotic solutions 
(where the top Lyaounov exponent Ai has positive values). For smaller fi 
chaotic solutions are grouped in the region of low values of an excitation 
frequency Q (Q < 0.27). Interestingly, the Lyapunov exponent, in that region 
Ai G [—0.03 : 0.03], is of the same ranges thus the change of fi does not scale 
them. Note this is a different region from that discussed by Chakraborty [12] 
who concentrated on the region where the resonance curve possessed multiple 
solutions. However in our particular system the chaotic solutions are evidently 
suppressed by higher excitation frequencies (Figs 3a-b). Thus continuing our 
research on chaotic solutions appearance in the Ueda system (Eq. 1) (Figs. 
2c-d and 3b) we will focus on low frequency region in this paper. 

Firstly we will look for homoclinic orbits which can be treated analytically 
by perturbation methods, namely by the Melnikov method. Such a treatment 
has been applied to selected problems in science and engineering [2,3]. 
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Fig. 4. Analytical amplitude C, calculated by harmonic balance approximation 
(Eq. 7), versus G [0.1 : 1.9] in the lowest order for system parameters: 7 = 1.0 
a = 0.05, and three values of /J, (/J, = 0.05, 0.10 and 0.15). 

2 Approximate solution around the main resonance 



In the vicinity of main resonance we assume periodic synchronized solution 

x = A{t)sm{Qt) + B(t)cos(ttt). (3) 

introducing it to Eq. 1 and making use of the following trigonometric identities: 



1 3 

(cos ip) 3 = — cos(3?/>) + - cos ip, 



1 3 

(sin-?/>) 3 = — - sin(3^>) + - sinip, 



sin("0) (cos if)) 2 = sin(?/>) — sin(?/>) 3 = - sin(3^>) + - sin if), 
cos(^)(sin^) 2 = cos("0) — cos(ip) 3 = — - cos(3^) + - cosip 



(4) 



we get 



d 2 A o2 , dA 
— — - Q 2 A+—a 
dt 2 dt 



aQB 



sm(Qt) 



7 (^ 3 -^5 2 )sin(3(^)) 



(5) 
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Fig. 5. Effective potential V e ff = V e ff(A,B) for the following system parameters: 
7 = 1.0 [i = 0.1, a = 0.05 and f2 = 0.21 plotted in two different scales (a) and (b) 
respectively. 



^| - Q 2 B + ^-a + aQA + Ab* + BA 2 ) 
dt 2 dt 4 



cos(fit) 



7(-^ 3 + ^^ 2 )cos(3(fit)). 



In the spirit of the harmonic balance approximation [10] we find fixed points 
neglecting higher harmonics sin(3f2t) and cos(3fit) in the lowest approxima- 
tion. Thus for A = and B = 



Of A - aVLB + -7 [A 6 + AB z j = //, 
-Q 2 B + aQA + ^7 (fi 3 + BA 2 ) = 0. 



(6) 



after some simple algebra we get simple equation 



c 2 a 2 n 2 + c 2 ( -n 2 + ^ 7 c 2 ) 2 - /i 2 = o, 



(7) 



where 



C = A + B 



The result of an analytical solution (Eq. 7) for the amplitude C versus fre- 
quency Q has been shown in Fig. 4. One can see the characteristic incline in 
the resonance in the right hand side. Above Q = 0.6 there are triple solutions 
where the upper ad bottom ones are stable and the middle one is unstable. 
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3 Melnikov approach beyond harmonic balance 



Motivated by Chakraborty [12] we now going beyond the harmonic balance 
approximation keeping all terms of (Eq. 4) Note, that the higher harmonic 
terms with 3Qt can be easily transformed into 2Qt and Qt respectively 

cos(3fit) = (2 cos(2fit) - 1) cos(ftt), (9) 
sin(3fit) = (2 cos(2fit) + 1) sin(fii). 

Thus introducing Eqs. 9 into Eq. 5 and simplifying by sin(f2t) and cos(f2t) we 
get 

^A-n 2 A + ^a-anB + 1 A 3 - f i = 1 ( -A 3 - -AB 2 ) cos(2fti), 
dt 1 dt \2 2 / 

, ] " B _ n 2 B + ^a + atlA + 1 B 3 = 7 (--B 3 - -BA 2 ) cos(2fit). 
dt V 2 2 / 

(10) 



dt 2 



In this way we have obtained new equations of motion for new coordinates 
A(t) and B(t) with a parametric excitations. 

Defining velocities an d vsit) and introducing small parameters into the 

equation e and corresponding parameters 7 and a (ej — 7, ea — a) the first 
order equations of motion as 



dVA 

dt 
dA 

~dt 



dv B 
dt 
dB 

~dt 



■■ Q 2 A - 7A 3 + e (-6X4 + aB + 7 {^A 3 - ^AB 2 ^j cos(2fit)) , 

■v A (11) 

= fi 2 5 - 7A 3 + e (-au B - aA + 7 (~^ 3 - ^A 2 B^J cos(2fit)) , 
= «b. (12) 



The effective unperturbed Hamiltonian for the above set of equations (Eq. 12) 
can be written 

H eff = y + y + V eff (A, B), (13) 
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Fig. 6. Homoclinic orbits in (B, Vb) plane for Vt = 0.8, 0.5, 0.2, respectively 
where 

Veff(A) = + ± T - ll A + — ^— + 1^, (14) 

and va,vb are effective velocities which define kinetic terms of the Hamiltonian 
Eq. 13. 

This potential for the chosen system parameters, 7 = 1.0 \i = 0.1, a = 0.05 
and Q = 0.21, has been plotted in Figs. 5a and b. Note that for more precise 
mesh (Fig. 5b) we observe double-well structure of potential with degenerated 
minima energy and a saddle point between them. 

Existence of this point with a horizontal tangent makes homoclinic bifurca- 
tions of the system possible i.e. transition from a regular to chaotic solution. 

Note the characteristic saddle point [A, B] = [A , 0) is going to be reached in 
exactly defined albeit infinite time t corresponding to +00 and —00 for stable 
and unstable orbits, respectively. On the other hand A can be obtained as 
the equilibrium fixed point from Eq. 11 and e = 

Q 2 , a , s 

A A - - = 0. 15 

7 7 

Using the Cardano formula 

^fF^ + fF^ (ie) 
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Fig. 7. Set of Melnikov curves plotted in a plane of system (u, 7). 'l'-'6' corresponds 
to fj, = 0.01, 0.03, 0.05, 0.08, 0.10 and 0.12, respectively. C and R symbolize regions 
denned by the critical Melnikov curve of possible chaotic and regular solutions. 

The homoclinic orbit can be derived by assuming that 

A = A and lim Bit) = 0. (17) 



In this case potential part for B variable reads 

V'b = -%B* + \b\ (18) 



Assuming the condition A = A and noting that Vb = at the saddle point 
(B = 0) we perform standard analysis on energy conservation formula 

0= V f + V B (B) (19) 



and after integration we get 

= t + t , (20) 
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Fig. 8. Phase portraits and Poincare maps for system parameters chosen as follows: 
ft = l.l, p = 7.0 (Fig. 8a) and /j = 7.5 (Fig. 8b); ft = 0.21, p = 0.01 (Fig. 8c) 
and fi = 0.35 (Fig. 8d). The initial condition used in both cases x(t = 0) = 0.0, 
v (t = 0) = 2.5. The top Lyapunov exponent Ai = 0.000, -0.025, -0.024, -0.024 for 
Figs 8a-d respectively. 



where to is a time like integration constant. 



Finally the homoclinic orbit is given by time dependent coordinate B(t) = 
B*(t) 

B * {t) = ^cosh(n(t + t )) (2i; 



and a corresponding velocity vs(t) = v* B {t) 
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Fig. 9. Top Lyapunov exponent against quadratic potential term 5 as in V S 2 (Eq. 2). 
For the smallest (negative) 5 the initial conditions we assumed to be x(t = 0) = 0.0, 
v(t = 0) = 2.5 while for any next larger 5 the final state values of x and v have been 
used as the initial conditions. 

The homoclinic orbits for the same system parameters 7 and fi as in Fig. 5 
and assumed three values of excitation frequency (Q = 0.8, 0.5, 0.2) has been 
plotted in Fig. 6. Interestingly, for higher Q we observe the effect of blowing 
the length of homoclinic loops and corresponding areas inside. 

In case of perturbed orbits IV s and W u the distance between them is given 
by the Melnikov function M(t ): 



+00 



M(t ) 



h(B*y B )Ag(B*,v* B )dt 



(23) 



where the corresponding differential forms h as the gradient of unperturbed 
Hamiltonian (Eq. 13) (for A(t) = A ) leading to equations of motion 

aF = -" B < T^ B =B ' (24) 

while g as its perturbation form of the above 



h = (n 2 B + ^B 6 )dB + v B dv B , (25) 
g = (-av B - aVLA + 7 (~^B 3 - ^BA 2 ^j cos(2fii)) dx 
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are defined on homoclinic manifold (B,vb) = (B*,v B ). Finally the Melnikov 
integral is given by 



+00 

M (t ) = J dt (-av* B - aQA + 7 (~^B* 3 - ^B*A* 2 ^j cos(2fit) 



(26) 



After substituting B*(t) and v* B (t) by formulae given in Eq. 6 and A*(t) = A 
(Eq. 16) we get (see Appendix A) 



4 Q 3 16 fl 4 
M (t ) = --a— + — 7— -Tre - * sin(2fit ) 
37 3 7^7 

Q 2 

+ \1*fA\ — ire' 7 " sin(2fit )- 

7 



(27) 



In the Melnikov theory M ~ d [1,2,3], where d is the distance between stable 
and unstable manifolds. Simple zero of the Melnikov function M is associ- 
ated with the cross-section of these manifolds indicating global homoclinic 
bifurcation. In our case this condition (for the set of parameters: a, 7, Q and 
A which a function of /i, 7 and A (/j,, 7, Q) (see Eqs. 15-16) is fulfilled for 
I sin(2Oi )| = 

an = (12Q 2 ^/7 + 9Aly/T)ne-* . (28) 



The analytical results for critical parameters 7 = 7 C and = fi c basing on 
this equation are shown in Fig. 7 for // = 0.01, 0.03, 0.05, 0.08, 0.10 and 0.12. 
respectively. Left and right handed regions denoted for each curve as C and 
R respectively are related to possible " Chaotic" and " Regular" solutions. The 
results show that for larger 7 chaotic solution is more limited. This tendency 
is better visible for a larger excitation amplitude [i. Note the formerly inves- 
tigated chaotic solution has been indicated by a singular point in the figure. 
Account for that case \x was chosen as 0.1 one can easily see that this solution 
match with the analytical prediction. We have also checked that our system 
undergoes typical period doubling cascade showing also three points solution 
as well (defined on Poincare maps) Figs. 8b, d. In case of ) This regulars solu- 
tion often accompany a chaotic solution. Note also that for one of presented 
solutions (in Fig. 8a) the Lyapunov exponent was approximately indicating 
a doubling period bifurcation point. A typical solution, synchronized with an 
excitation term, has been shown in Fig. 8c. 
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4 Summary and Conclusions 

In summary we have performed the Mielnikov analysis for the Ueda system 
with a single well potential. Through transforming the system to new variables 
it was possible to investigate the Mielnikov criterion for a global homoclinic 
bifurcation from regular to chaotic oscillations. 

Our investigation was limited purposely to a small frequency Q where the 
chaotic solutions emerge numerically (Figs. 2c, d). However for different non- 
liner systems involving nonlinear damping terms and the self-excitation effects 
[13] the region of chaotic solutions could be different. The main simplification 
in our treatment was the assumption A(t) = A =const. In a more general 
case one should expect additional time dependence of the amplitude A which 
could create an additional shift of the critical lines in Fig. 7. This shift should 
be dependent on Q, which influences strongly the size of homoclinic orbits 
Fig. 6. Consequently, for our system, in the limit of large Q the size of the 
homoclinc orbit is so large that condition A(t) =const. cannot be applied. 

Our results for the top Lyapunov exponent (Fig. 9) show also that the chaotic 
solution was preserved in presence of a small linear force term 5 (see V S 2 in 
Eq. 2) so the method presented here can be generalized to the Duffing system 
having linear and cubic force terms. 
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Appendix A 

After substituting B*(t), v* B (t) by formulae given in Eq. 6 and assuming that 
A*(t) = A as calculated in Eq. 16 into Eq. 23, taking r = \/5t/2 we get 





(A.l) 
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+ 3 , A 2^_ T sinh(^(t + t )) (2mdt 



7 J cosh 6 (Q(t + t )) 
Using r = Q(t + to) we obtain 



+00 9 /— +00 

, ,/ x ~2^o /■ sinh t , , ( V2„, f smhr , 

7 7 cosh r x /7 7 cosh r 



2 4 f sinhr 



+ 7 Vl A / — -^cos(2r-2fit )dt (A.2) 



7 \fl 7 cosh r 
v —00 



3 /* smh 7~ 
+ ^A -n 2 / — -^cos(2r-2fit )dt. 
7 7 cosh r 



7 7 cosh r 

—00 

The above Melnikov integral (Eq. A.2) can be written as 



2 o3r 1 ~ A ^n2 



M(t ) = -a-tt 6 h + aA—tt 2 I 2 

7 ^fl 

+ 7-^=^3 + lA-Q 2 h. (A.3) 

7v /7 7 



Integrals ii and I 2 can be calculated directly 

h = | ^2 = 0. (A.4) 

Let us write I3 and I4 in the complex space as 



+00 . +00 . 

/ ^cos(2r)dt + i / e— sin(2r)dn 

7 cosh t J cosh r J 

—00 —00 

-(-oo "4~oo 

+ sin(2^t )Im| J cos(2r)dt + z / ~^5~ sin(2r)dt|, (A.5) 



and 



hoo . +00 . 

{/" sinh t f sinh t 1 
/ o— cos(2r)dt + i / o— sin(2r)cit 
7 cosh r 7 cosh r J 
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f (' Slllll T f Slllll 7~ 1 

+ sin(2fit )Im<^ / cos(2r)dt + i / =- sin(2r)rft L (A.6) 

I J cosh r j cosh t J 



-oo — oo 



respectively. 

Now we can simplify the notation 



J 3 = cos(2fit )Re( / ^-^e 2iT dt\ 
I J cosh T J 



-oo 



+ sin(2^ )Im<| / ^-^e 2lT dt \, (A.7) 
' cosh r 

sinh' 



-oo 
hoc 



{[' sinh t 1 
/ ^e 2lT dt\ 
J cosh r J 

— oo 
+00 . 

f c sinn t 

+ sin(2fit )Im / ^e 2lT dt \. (A.8) 

I j cosh r J 



Applying the residue theorem, we get 



Is = cos(2f2t )Re ^-7n(cosh -n — sinh 7r)^ 



+ sin(2nt )Im ^-7ri(cosh n — sinh rr) j , (A. 9) 

^4 = cos(2f2t )R.e(4"7ri(cosh 7r — sinh 7r)) 

+ sin(2^t )Ini(47ri(cosh7r - sinliTr)). (A.10) 



Consequently 

g 

J 3 = sin(2f2t )-7r(cosh7r — sinhvr), (A. 11) 

o 

J 4 = sin(2fit )47r(cosh7r - sinli7r). (A. 12) 

Finally the Melnikov integral reads 



4 fi 3 16 fi 4 
M (t ) = + — 7^-7re- n sin(2fit ) 

37 3 7^ 

+ I27A 2 — vre"^ sin(2fit )- (A.13) 

7 
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